# NIPT Y染色体浓度分析问题一完整解决方案
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

class NIPTAnalyzer:
    def __init__(self):
        self.boy_data = None
        self.girl_data = None
        self.analysis_results = {}
        
    def load_data(self):
        """加载NIPT数据"""
        print("正在加载NIPT数据...")
        try:
            self.boy_data = pd.read_excel('data/附件.xlsx', sheet_name='男胎检测数据')
            self.girl_data = pd.read_excel('data/附件.xlsx', sheet_name='女胎检测数据')
            print(f"成功加载数据: 男胎{len(self.boy_data)}条, 女胎{len(self.girl_data)}条")
            return True
        except Exception as e:
            print(f"数据加载失败: {e}")
            return False
    
    def explore_data(self):
        """数据探索"""
        print("\n=== 数据探索 ===")
        
        # 查看男胎数据结构
        print("男胎数据列名:")
        boy_cols = self.boy_data.columns.tolist()
        for i, col in enumerate(boy_cols):
            print(f"{i+1}. {col}")
        
        # 关键字段检查
        key_fields = ['Y染色体浓度', '孕妇BMI']
        available_fields = [f for f in key_fields if f in boy_cols]
        missing_fields = [f for f in key_fields if f not in boy_cols]
        
        print(f"\n可用字段: {available_fields}")
        if missing_fields:
            print(f"缺失字段: {missing_fields}")
        
        return available_fields
    
    def analyze_y_bmi_relationship(self):
        """分析Y染色体浓度与BMI关系"""
        print("\n=== Y染色体浓度与BMI关系分析 ===")
        
        # 数据清理
        required_cols = ['Y染色体浓度', '孕妇BMI']
        boy_clean = self.boy_data[required_cols].dropna()
        
        print(f"有效数据点数: {len(boy_clean)}")
        
        # 基本统计
        print("\nY染色体浓度统计:")
        print(boy_clean['Y染色体浓度'].describe())
        
        print("\n孕妇BMI统计:")
        print(boy_clean['孕妇BMI'].describe())
        
        # 相关性分析
        correlation, p_value = stats.pearsonr(boy_clean['Y染色体浓度'], 
                                           boy_clean['孕妇BMI'])
        
        print(f"\nPearson相关系数: {correlation:.4f}")
        print(f"P值: {p_value:.4f}")
        
        # 显著性判断
        if p_value < 0.001:
            significance = "极显著"
        elif p_value < 0.01:
            significance = "非常显著"
        elif p_value < 0.05:
            significance = "显著"
        else:
            significance = "不显著"
        
        print(f"相关性: {significance}")
        
        return boy_clean, correlation, p_value
    
    def build_regression_model(self, data):
        """建立回归模型"""
        print("\n=== 回归模型建立 ===")
        
        X = data['孕妇BMI'].values.reshape(-1, 1)
        y = data['Y染色体浓度'].values
        
        # 线性回归
        linear_model = LinearRegression()
        linear_model.fit(X, y)
        
        # 预测和评估
        y_pred = linear_model.predict(X)
        r2 = r2_score(y, y_pred)
        
        # 模型显著性检验
        n = len(data)
        k = 1  # 自变量个数
        df_total = n - 1
        df_residual = n - k - 1
        
        # F检验
        ss_total = np.sum((y - np.mean(y))**2)
        ss_residual = np.sum((y - y_pred)**2)
        ss_regression = ss_total - ss_residual
        
        ms_regression = ss_regression / k
        ms_residual = ss_residual / df_residual
        f_stat = ms_regression / ms_residual
        
        # 计算p值
        p_value_f = 1 - stats.f.cdf(f_stat, k, df_residual)
        
        print(f"回归方程: Y = {linear_model.intercept_:.4f} + {linear_model.coef_[0]:.4f} * BMI")
        print(f"R² = {r2:.4f}")
        print(f"F统计量 = {f_stat:.4f}")
        print(f"P值 = {p_value_f:.4f}")
        
        return {
            'model': linear_model,
            'r2': r2,
            'f_stat': f_stat,
            'p_value': p_value_f,
            'equation': f"Y = {linear_model.intercept_:.4f} + {linear_model.coef_[0]:.4f} * BMI"
        }
    
    def bmi_group_analysis(self, data):
        """BMI分组分析"""
        print("\n=== BMI分组分析 ===")
        
        # 按BMI分组
        bins = [0, 20, 28, 32, 36, 40, 100]
        labels = ['<20', '[20,28)', '[28,32)', '[32,36)', '[36,40)', '≥40']
        data['BMI_group'] = pd.cut(data['孕妇BMI'], bins=bins, labels=labels, right=False)
        
        # 各组统计
        group_stats = data.groupby('BMI_group')['Y染色体浓度'].agg([
            'count', 'mean', 'std', 'min', 'max'
        ])
        
        print("\n各BMI组Y染色体浓度统计:")
        print(group_stats)
        
        # 方差分析
        groups = [group['Y染色体浓度'].values for name, group in data.groupby('BMI_group')]
        f_stat, p_value = stats.f_oneway(*groups)
        
        print(f"\n方差分析:")
        print(f"F统计量: {f_stat:.4f}")
        print(f"P值: {p_value:.4f}")
        
        return group_stats, f_stat, p_value
    
    def create_analysis_report(self):
        """创建分析报告"""
        print("\n=== 生成分析报告 ===")
        
        # 加载和探索数据
        if not self.load_data():
            return
        
        available_fields = self.explore_data()
        
        if len(available_fields) < 2:
            print("数据缺少必要字段，无法完成分析")
            return
        
        # 主要分析
        clean_data, corr, p_corr = self.analyze_y_bmi_relationship()
        regression = self.build_regression_model(clean_data)
        group_stats, f_stat, p_anova = self.bmi_group_analysis(clean_data)
        
        # 创建可视化
        self.create_visualizations(clean_data, regression)
        
        # 总结报告
        self.generate_summary_report(clean_data, regression, group_stats)
    
    def create_visualizations(self, data, regression):
        """创建可视化图表"""
        fig, axes = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle('NIPT Y染色体浓度与BMI关系分析', fontsize=16)
        
        # 1. 散点图与回归线
        ax1 = axes[0, 0]
        ax1.scatter(data['孕妇BMI'], data['Y染色体浓度'], alpha=0.6, s=20)
        ax1.set_xlabel('孕妇BMI')
        ax1.set_ylabel('Y染色体浓度')
        ax1.set_title('散点图与回归线')
        
        # 添加回归线
        x_line = np.linspace(data['孕妇BMI'].min(), data['孕妇BMI'].max(), 100)
        y_line = regression['model'].predict(x_line.reshape(-1, 1))
        ax1.plot(x_line, y_line, 'r-', linewidth=2, label=f'R²={regression["r2"]:.3f}')
        ax1.legend()
        
        # 2. 残差图
        ax2 = axes[0, 1]
        X = data['孕妇BMI'].values.reshape(-1, 1)
        y = data['Y染色体浓度'].values
        y_pred = regression['model'].predict(X)
        residuals = y - y_pred
        
        ax2.scatter(y_pred, residuals, alpha=0.6, s=20)
        ax2.axhline(y=0, color='r', linestyle='--')
        ax2.set_xlabel('预测值')
        ax2.set_ylabel('残差')
        ax2.set_title('残差分析图')
        
        # 3. BMI分组箱线图
        ax3 = axes[1, 0]
        data.boxplot(column='Y染色体浓度', by='BMI_group', ax=ax3)
        ax3.set_xlabel('BMI分组')
        ax3.set_ylabel('Y染色体浓度')
        ax3.set_title('各BMI组Y染色体浓度分布')
        
        # 4. 直方图
        ax4 = axes[1, 1]
        ax4.hist(data['Y染色体浓度'], bins=30, edgecolor='black', alpha=0.7)
        ax4.set_xlabel('Y染色体浓度')
        ax4.set_ylabel('频数')
        ax4.set_title('Y染色体浓度分布')
        
        plt.tight_layout()
        plt.savefig('NIPT_Y染色体浓度分析.png', dpi=300, bbox_inches='tight')
        plt.show()
    
    def generate_summary_report(self, data, regression, group_stats):
        """生成总结报告"""
        print("\n" + "="*50)
        print("NIPT Y染色体浓度分析问题一分析报告")
        print("="*50)
        
        print("\n【数据概况】")
        print(f"- 总样本数: {len(data)}")
        print(f"- BMI范围: {data['孕妇BMI'].min():.1f} - {data['孕妇BMI'].max():.1f}")
        print(f"- Y染色体浓度范围: {data['Y染色体浓度'].min():.3f} - {data['Y染色体浓度'].max():.3f}")
        print(f"- Y染色体浓度均值: {data['Y染色体浓度'].mean():.3f}")
        
        print("\n【相关性分析】")
        correlation, p_value = stats.pearsonr(data['Y染色体浓度'], data['孕妇BMI'])
        print(f"- Pearson相关系数: {correlation:.4f}")
        print(f"- P值: {p_value:.4f}")
        
        if p_value < 0.05:
            print("- 结论: Y染色体浓度与BMI存在显著相关性")
        else:
            print("- 结论: Y染色体浓度与BMI相关性不显著")
        
        print("\n【回归模型】")
        print(f"- 模型类型: 线性回归")
        print(f"- 回归方程: {regression['equation']}")
        print(f"- 决定系数R²: {regression['r2']:.4f}")
        print(f"- F统计量: {regression['f_stat']:.4f}")
        print(f"- 模型显著性P值: {regression['p_value']:.4f}")
        
        if regression['p_value'] < 0.05:
            print("- 结论: 回归模型在统计上显著")
        else:
            print("- 结论: 回归模型在统计上不显著")
        
        print("\n【BMI分组分析】")
        print("各组Y染色体浓度统计:")
        print(group_stats)
        
        print("\n【临床意义】")
        print("- 根据4%的Y染色体浓度阈值标准:")
        high_conc = data[data['Y染色体浓度'] >= 0.04]
        low_conc = data[data['Y染色体浓度'] < 0.04]
        print(f"  高浓度(≥4%): {len(high_conc)}例 ({len(high_conc)/len(data)*100:.1f}%)")
        print(f"  低浓度(<4%): {len(low_conc)}例 ({len(low_conc)/len(data)*100:.1f}%)")
        
        print("\n【建议】")
        print("- 基于分析结果，建议根据孕妇BMI值制定个性化的NIPT检测策略")
        print("- 高BMI孕妇可能需要调整检测时点以提高准确性")

# 主程序
if __name__ == "__main__":
    analyzer = NIPTAnalyzer()
    analyzer.create_analysis_report()